rm(list = ls())
library(maps)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(mapdata)
library(ggplot2)
source("../R/geomatching.R")
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
source("../R/idw2hr.R")
## Loading required package: iterators
## Loading required package: parallel
source("../R/hr2poly.R")
AQ_EEA_NO2: NO2 point data measured from air quality monitoring data.
AQ_CAMS_N02: NO2 gridded data (0.01° x 0.01°) obtained from CAMS mathematical model.
# CSV data (.csv)
AQ_EEA_NO2 <- read.csv("AQ_EEA_NO2.csv")
# R environment (.RData)
load("AQ_CAMS_NO2.RData")
geomatchingdata_dict <- list(
"AQ_EEA_NO2" = list(
"data" = AQ_EEA_NO2,
"format" = "xyt",
"type" = "points",
"crs" = 4979
),
"AQ_CAMS_NO2" = list(
"data" = AQ_CAMS_NO2,
"format" = "matrix",
"type" = "grid",
"crs" = 4326
)
)
data_ls <- list()
settings = list(
"format"=c(),
"type"=c(),
"crs"=c()
)
for (i in seq_along(data_dict)) {
data_ls[[i]] <- data_dict[[i]]$data
settings$format[[i]] <- data_dict[[i]]$format
settings$type[[i]] <- data_dict[[i]]$type
settings$crs[[i]] <- data_dict[[i]]$crs
}
res_geomatch <- geomatching(
data=data_ls,
settings = list(
"format"=settings$format,
"type"=settings$type,
"crs"=settings$crs
)
)
## [1] "reading data 1"
## [1] "done"
## [1] "reading data 2"
## [1] "done"
## [1] "converting data 1 to ST"
## [1] "done"
## [1] "converting data 2 to ST"
## [1] "done"
# rename columns
res_geomatch <- res_geomatch %>%
rename(
altitude = Altitude,
AQ_CAMS_NO2 = matrix_2
)
for (var in c("AQ_EEA_NO2", "AQ_CAMS_NO2", "altitude")) {
p <- ggplot() +
geom_map(
data = map_data("italy"),
map = map_data("italy"),
aes(map_id = region),
fill = "lightgrey",
color = "black",
linewidth = .1,
alpha=.9
) +
geom_point(
data = res_geomatch[res_geomatch$time==res_geomatch$time[1],],
aes(x = longitude, y = latitude, color = .data[[var]]),
size = 5,
) +
labs(
title = sprintf("Geographical map of %s", var),
subtitle = "Output of geomatching",
x = "Longitude",
y = "Latitude"
) +
scale_color_viridis_c(
option = "C"
) +
theme(
legend.title = element_blank(),
plot.title = element_text(face = "bold")
)
print(p)
}
idw2hrres_idw2hr <- idw2hr(
data = res_geomatch,
col_names = list(
"lon" = "longitude",
"lat" = "latitude",
"t" = "time"
),
interest_vars = list("altitude", "AQ_EEA_NO2", "AQ_CAMS_NO2"),
crs = 4979
)
for (var in c("AQ_EEA_NO2", "AQ_CAMS_NO2", "altitude")) {
p <- ggplot() +
geom_point(
data = res_idw2hr[res_idw2hr$time==res_idw2hr$time[1],],
aes(x = longitude, y = latitude, color = .data[[var]]),
size = 5,
) +
geom_map(
data = map_data("italy"),
map = map_data("italy"),
aes(map_id = region),
fill = "lightgrey",
color = "black",
linewidth = .1,
alpha=.05
) +
labs(
title = sprintf("Geographical map of %s", var),
subtitle = "Output of idw2hr",
x = "Longitude",
y = "Latitude"
) +
scale_color_viridis_c(
option = "C"
) +
theme(
legend.title = element_blank(),
plot.title = element_text(face = "bold")
)
print(p)
}
hr2polyres_hr2poly <- hr2poly(
res_idw2hr,
polygon_type = "mun",
crs = 4979
)
for (var in c("AQ_EEA_NO2_min", "AQ_EEA_NO2_mean", "AQ_EEA_NO2_median",
"AQ_EEA_NO2_max", "AQ_EEA_NO2_std", "AQ_CAMS_NO2_min", "AQ_CAMS_NO2_mean",
"AQ_CAMS_NO2_median", "AQ_CAMS_NO2_max", "AQ_CAMS_NO2_std", "altitude_min",
"altitude_mean", "altitude_median", "altitude_max", "altitude_std")) {
p <- ggplot() +
geom_sf(
data=res_hr2poly,
aes_string(fill=var)
) +
labs(
title = sprintf("Geographical map of %s", var),
subtitle = "Output of hr2poly",
x = "Longitude",
y = "Latitude"
) +
scale_fill_viridis_c(
option = "C"
) +
theme(
legend.title = element_blank(),
plot.title = element_text(face = "bold")
)
print(p)
}
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.